Self-consistent Skyrme QRPA for use in axially-symmetric nuclei of arbitrary mass 
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We describe a new implementation of the quasiparticle random phase approximation (QRPA) 
in axially-symmetric deformed nuclei with Skyrme and volume-pairing energy-density functionals. 
After using a variety of tests to demonstrate the accuracy of the code in ^^'^^Mg and ^^O, we report 
the first fully self-consistent application of the Skyrme QRPA to a heavy deformed nucleus, calcu- 
lating strength distributions for several in ^^^Yb. We present energy-weighted sums, properties 
of 7-vibrational and low-energy states, and the complete isovector El strength function. 

The QRPA calculation reproduces the properties of the low-lying 2+ states as well or better than 
it typically does in spherical nuclei. 
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I. INTRODUCTION 

The quasiparticle random phase approximation 
(QRPA) ^ has a long history in nuclear physics. Its 
virtues include applicability to many types of excitation 
across the isotopic chart, preservation of energy- weighted 
sum rules, and elimination of spurious motion. In addi- 
tion, the QRPA has several appealing interpretations; it 
is both a boson approximation for collective modes and 
the small-amplitude limit of the time-dependent Hartree- 
Fock-Bogoliubov (HFB) approximation. Its downside, 
traditionally, has been a limited ability to describe large- 
amplitude motion and complicated non-collective states, 
deficiencies that prompted the development of several 
more complicated methods, as, e.g., in Refs. 0, 

Recent years, however, have seen a revival of the 
QRPA, despite its drawbacks. The primary reason is the 
increasing connection between nuclear mean-field theory 
and density-functional theory (DFT) [H, Q . The notion 
that Hartree-Fock or HFB calculations can be relevant 
beyond their naive range of validity has motivated at- 
tempts to describe a wide-range of nuclear properties 
in mean-field theory and extensions. The QRPA is the 
most straightforward extension that fits into the DFT 
paradigm; to the extent that the energy functional used 
in HFB calculations is exact, the QRPA provides the ex- 
act linear (i.e. small-amplitude) response function in the 
adiabatic limit Combined with its other features, 
its connection with DFT makes the QRPA an important 
tool in attempts both to develop a "Universal Nuclear 
Energy-Density Functional" (UNEDF) , and to apply the 
functional to, e.g. nuclear astrophysics Q. 

The prototype energy-density functional is of Skyrme 
form, corresponding roughly to effective interactions that 
have zero range, with derivatives simulating finite-range 
effects. In the last five or ten year, a number of groups 
have developed self-consistent (Q)RPA codes for use with 
these functionals or similar finite-range and relativistic 
versions, first in spherical nuclei, (see Ref. 3 and refer- 
ences therein), and then in axially symmetric deformed 
nuclei [ol-fTsj. Heavy deformed nuclei are still problem- 
atic, however. Though the deformed RPA, without pair- 



ing, is now tractable in heavy nuclei jl4l |. a separable 
approximation to the Skyrme- QRPA equations has been 
applied in such nuclei [isl IT6j . and efficient new meth- 
ods for solving the full QRPA equations are promising 
[3, [l3| , the numerical complexity of deformed systems 
has so far limited fully self-consistent Skyrme-QRPA cal- 
culations to nuclei with A < AO. In this paper, we present 
a highly parallelized version of the Skyrme QRPA that 
we are beginning to apply to heavy deformed nuclei. Af- 
ter discussing the structure of the code and demonstrat- 
ing its accuracy, we present a preliminary application to 
the nucleus ^^^Yb. We focus on results here, postponing 
most of the formalism to a later publication. 

Our QRPA code comes in several versions, developed 
successively as we progressed from tests in light nuclei 
to fuU-fiedged calculations in heavy nuclei. All versions 
treat the entire Skyrme -|- Coulomb functional in a way 
that is completely consistent with HFB calculations (re- 
stricted for the time being to "volume" pairing). In ad- 
dition, all the versions preserve axial and parity sym- 
metries and the time-reversal invariance of the ground 
state, and therefore require an HFB code that produces 
single-quasiparticle wave functions of the two variables 
r = + and z, with M = {j^) a good quan- 

tum number. Finally, all diagonalize the traditional 
QRPA A, B matrix [l| in the "M-scheme", and use 
the rigid- rotor approximation [l^ [3l i with the deformed 
QRPA solutions as intrinsic states, to calculate transi- 
tion strength. The versions differ, however, in the ba- 
sis in which the HFB equations are solved, in the basis 
in which the QRPA matrix is constructed, and in the 
way wave functions are represented numerically. Ver- 
sion a) takes quasiparticle-basis wave functions from the 
transformed-oscillator code HFBTHO [l^ (though we 
use the ordinary harmonic-oscillator basis), represents 
the wave functions on an equidistant mesh, and con- 
structs the QRPA matrix in the quasiparticle representa- 
tion. Version b) substitutes the Vanderbilt "cylindrical- 
box" S-spline-based HFB code [U for HFBTHO. Ver- 
sion c) modifies the QRPA part of version b) by using the 
canonical-quasiparticle basis, represented with i?-splines, 
in place of the quasiparticle basis to speed calculation 
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and save memory. The vast majority of the computing 
time in all versions is in the construction of the QRPA 
matrix, each element of which requires a series of two- 
dimensional integrals for the Skyrme interaction, and an 
additional multipole expansion for the Coulomb interac- 
tion. The set of matrix elements can be divided among 
many thousands of processors so that the calculation is 
manageable on fast supercomputers. 

Section II below describes tests in relatively light sys- 
tems. Section III presents an application to heavy nuclei. 



II. TESTS 

To display the accuracy of our codes, we show the re- 
sults of several tests in nuclei with ^ < 40. We start 
with spurious states. A fully self-consistent and numer- 
ically perfect QRPA will completely separate spurious 
states from physical ones, and put them at zero energy. 
Small numerical errors can spoil the treatment of spu- 
rious states, however, so any calculation that separates 
them well has passed a serious test. 

Figure [T] shows transition matrix elements of the num- 
ber operators, 



5,- = |(0|iV,|fc)|^ 



(1) 



with k labeling excitations and r protons/neutrons, to 
states with K"^ = 0+ in a calculation of ^^Mg. To ob- 
tain the QRPA wave functions, we use version a) of our 
code with the Skyrme interaction SkP j22| . restricting 
ourselves to three harmonic oscillator shells. The HFB 
calculation yields pairing gaps of Ap = 1.681 MeV (pro- 
tons), A„ = 1.426 MeV (neutrons) and no quadrupole 
deformation in this small single-particle space. This re- 
sult docs not impl y t hat a calculation in a larger space 
also gives (3 = [2j]. We use a small space to allow 
a complete QRPA treatment of two-quasiparticle excita- 
tions; a full separation of the spurious = 0"*" strength 
associated with the particle-number violation requires no 
less. And we indeed achieve essentially perfect separa- 
tion. The figure shows negligible strength to all excita- 
tions except the two spurious states, which come out at 
E = 0.008 McV and 0.024 MeV. 

Next, we turn to spurious rotation, arising in the K'^ 
= 1+ channel. Here the nucleus must be deformed. Our 
HFB solution, again with SkP but for ^^Mg in five oscil- 
lator shells, yields (3 = 0.28, Ap = 0.034 MeV, and A„ 
= 0.131 MeV. Figure [2] displays the resulting transition 
strengths for the operator J-i = — iJy Most of the 
strengths arc four order of magnitude smaller than that 
of the spurious state at E' = 0.045 McV. 

Finally, we examine spurious translational motion. 
Here, an oscillator basis for HFB is not optimal, and 
so we employ versions b) and c) of our code (they give 
identical strength functions) , which use wave functions in 
a large cylindrical box. For these calculations, in ^^Mg, 
we choose the box to have size r^ax = -^max = 10 fm 
and use the Skyrme functional SLy4 [2^. We take a 
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FIG. 1: (Color online) Transition strength for the proton and 
neutron number operators to = 0"'' states in ^^Mg with 
the Skyrme functional SKP (see text). 
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FIG. 2: (Color online) Transition strength for the operator 
J_i to A"" = 1+ states in ^''Mg with SkP. 
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FIG. 3: (Color onUne) IS El transition strength to if" = 0" 
states in ^®Mg with SLy4. Solid lines represent strength that 
is corrected to eliminate residual spurious contributions. 



3 



TABLE I: Comparison of energies and (corrected) IS El 
strengths for the three lowest-lying J'^ = 1~ states in ^®0, 
calculated with spherical (J-scheme) code and the current 
(M-scheme) code and SLy4. The correction barely changes 
the strength for the two physical excited states, but reduces 
that of the spurious (lowest) state by many orders of magni- 
tude. 



J-scheme, — 1 


M-scheme, K'' = 0" 


E 


oISBl 
^cor 


E 


'-'cor 


(MeV) 


(fm«) 


(MeV) 


(fm'') 


0.323 


7.051x10"^ 


0.472 


1.298x10^"* 


7.500 


1.461x10 


7.440 


1.433x10 


10.610 


5.739x10"^ 


10.681 


4.283x10"^ 



pairing strength V — —140.348 MeV fm'^ and include 
quasiparticle states with energy < 300 MeV in the den- 
sity and pairing tensor. These parameters yield /? = 
-0.27, Ap = 1.365 MeV, and A„ = 0.002 MeV. Figure 
[3] shows isoscalar (IS) E\ transition strengths, for which 
the excitation operator is X)iLi to states with 
K'^ = 0~. The figure contains two sets of lines, the sec- 
ond of which adds a correction term to the IS operator 
(via the prescription of Ref. [l^ ) to remove residual spu- 
rious strength from physical excitations. The difference 
between the two sets is very small in all the physical 
states shown, and completely negligible in higher-energy 
states. In sum, our code is accurate enough to handle 
spurious states in this mass region extremely well. 

Since we have done extensive calculations with a spher- 
ical J-scheme code over the past few years [1^, we can 
test our current codes further by comparing their results 
with those obtained in the J scheme. Table |T] shows 
energies and IS E\ transition strengths (with the cor- 
rection mentioned above) for the three lowest J'^ = 1^ 
levels, along with the corresponding K"^ = 0^ energies 
and strengths from version b) of our current code, in the 
spherical nucleus ^®0. The two codes take wave func- 
tions from entirely different HFB codes: a slightly modi- 
fied version of HFBRAD [H called HFBMARIO for the 
spherical QRPA and the Vanderbilt HFB code for the 
deformed QRPA. The first state on each side of the table 
is the spurious state, with very small strength because of 
the correction. The next two, both genuine excitations, 
are nearly the same in both energy and strength in the 
two calculations. The full strength function, folded with 
a Lorentzian of width 3 MeV (see in Eq. (1) in Ref. Q) 
displayed in Fig. |4l shows the same level of agreement. 
The very small differences in the continuum are due to 
differing box boundary conditions: the spherical calcu- 
lation is in a spherical box with radius 20 fm and the 
deformed calculation is in the same cylindrical box we 
used for ^^Mg. 

Before moving to heavy nuclei, we display the results 
of one more test — this time simultaneously checking 
the Vanderbilt HFB code [2l| underlying our deformed 
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FIG. 4: (Color online) Full IS E\ strength function corre- 
sponding to Tab. H] The letters J and (M) denotes the J- 
(AJ-) scheme calculations. 



QRPA, and our procedure for constructing the canonical 
basis by diagonalizing the density matrix (which is the 
same as in Ref. Q). Figure [5] displays the proton 0p3/2 
canonical basis wave function in ^^O produced by both 
the spherical and deformed procedures. The agreement 
is perfect and far from trivial. Though Ap = in this 
calculation, we still construct the density matrix from 
the Vanderbilt-HFB output and diagonalize it to obtain 
the deformed canonical wave function. (There is no arbi- 
trariness in this wave function because only one proton 
P3/2 state is occupied in Oxygen.) The other bound-state 
wave functions produced by the two procedures, though 
we don't display them, agree equally well. 
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FIG. 5: (Color online) Comparison of proton OP3/2 canonical- 
basis wave function produced by J- and M-scheme codes for 
^^O at 2 = 0. The labels s-dn and s-up denote the spin-down 
and spin-up components in the M-scheme calculation. In the 
J scheme, the spin-down component is identically zero. 
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TABLE II: Energy-weighted sums for strength functions, in 
our QRPA calculations and from analytical sum-rule expres- 
sions. We include QRPA states with up to 90 MeV of ex- 
citation energy. The units of the IS El sum are MeV fm®, 
and those of the isovector (IV) El sum are MeV fm^ Those 
of all E2 sums are MeV fm*. The IS El strength has been 
corrected to remove spurious components. The contribution 
from negative values of K is not included. 



Transition 




QRPA 


Analytical 


operator 


of solution 






IS El 


1" 


1043560 


1042413 


IV El 


1" 


289.819 


285.764 


IS El 


0" 


2015266 


2019465 


IV El 


0" 


291.859 


285.764 


IS E2 


2+ 


64700 


63877 


IV E2 


2+ 


20284 


20076 


IS E2 


1+ 


76159 


88197 


IV E2 


1+ 


28517 


28174 


IS E2 


0+ 


97886 


96867 


IV E2 


0+ 


31271 


30874 



III. HEAVY NUCLEI 

Having thoroughly tested several versions of the de- 
formed QRPA, we apply version c), representing the best 
combination of speed and accuracy, to the nucleus ^^^Yb. 
We set up the HFB calculation as follows: we use a 
"box" with Tmax = Zmax = 20 fm, cut off the quasi- 
particle spectrum at 60 MeV and take the maximum 
z-component of quasiparticle angular momentum to be 
19/2; these parameters define a single-quasiparticle space 
with 4648 proton states and 5348 neutron states. We use 
the Skyrme functional SkM* [26j with volume pairing 
strengths Vp = -218.521 MeV fm^ and K = -176.364 
MeV fm'^ (determined from measured odd-even mass dif- 
ferences). The calculation yields Ap = 1.248 MeV, A„ 
= 0.773 MeV, and /3 = 0.34. 

We input half the canonical wave functions, those 
with jz > 0, in the QRPA, and construct and truncate 
two-canonical-quasiparticle configurations in the follow- 
ing way: We define critical particle-particle and particle- 
hole occupation probabilities (^cut)^ = 10~^ and (Wcut)^ 
= 10~^". If both canonical states have occupation prob- 



abilities vf and w| such that vf,v^>l- (Wcut)^ 



or vt 



vj < (Wcut)^ 1 '^^ omit the configuration. We also omit 
configurations with vf < w| for which vf/v'j < (v^^^)'^. 
These cuts result in a QRPA matrix whose size, while de- 
pending on multipolarity, is typically about 160,000 by 
160,000. 

Tabic im shows energy- weighted sums, alongside values 
obtained from sum rules, for IS and isovector (IV) electric 
operators in all ^^^Yb channels that we calculate. The 
differences between the QRPA sums and the sum rules 
are less than about 2% except in the IS K'^ = 1+ channel. 



TABLE III: Energies and B{E2;0+ 2+)'s for the 7- 
vibrational and "/?- vib" states of ^^^Yb. Experimental data 
are from Ref. [2^ (see also [1^). For the definition of B{E2), 
see, e.g., Ref. p^ |. 







Exp. 


Cal. 


7- vib. 


E (MeV) 


1.466 


2.261 




B{E2) (e^b^) 


0.0433 -f-5-15 


0.041 


"/3-vib." 


E (MeV) 


1.117 


1.390 




B{E2) (e^b^) 


0.0081 17 


0.00495 



where we were unable to adequately separate the spurious 
rotational state without going to a larger space. In other 
channels the spurious mode is under better control. In 
the IS 0^ and 1" channels, we have some contamination 
at low energies, but it is weak enough that the subtrac- 
tion procedure of Ref. [l2| restores the sum rule nearly 
exactly. In the 0"^ channel, as the tabic shows, the sep- 
aration is quite good even without subtraction. The 1^ 
channel can be corrected as well, but doing so requires a 
numerical procedure that we have not yet implemented. 

The accuracy of the energy- weighted sums in the K'^ = 
0+ and 2+ channels indicates that our approach is reli- 
able for low-lying quadrupole shape vibrations. Table 
Hm shows the energies and B{E2;0+ 2+)'s of both 
the 7-vibrational K'^ — 2+ state and a low-energy I-C^ 
= 0+ state with a significant B{E2) that we denote by 
"/3-vib". The quotation marks indicate that, unlike the 
clearcut 7-vibrational state, the 0"*" state has a somewhat 
smaller B{E2) than is typical of vibrational modes. Both 
states have been studied experimentally, e.g. in Ref. [27| . 

The agreement of both the energies of these states and 
their transition strengths with measured values is at a 
level that is typical of QRPA calculations in spherical 
nuclei. In Ref. [s^l we investigated a large set of such 
nuclei, characterizing the quality of the QRPA by two 
quantities: 



Re = HEc.i/E, 



cxp; ) 



(2) 



Rq = In JB{E2),a,/B{E2) 



cxp J 



where sufiices cal and exp denote calculated and exper- 
imental. The results in Tab. IIIII correspond to Re = 
0.43 and Rq = -0.03 for the 7-vibration, and Re = 0.22 
and Rq = —0.25 for "/3-vib." The histograms in Figs. 4 
and 9 of Ref. [s^ show these values to be near the most 
common values in spherical nuclei. 

Figure [H] displays the IV El strength function. The 
thick curve is the sum of strengths in all channels, and 
can be compared with experimental data. The peaks of 
the K'^ = 0~ and 1~ distributions in Fig.[6]lie at different 
energies, as is often the case in deformed nuclei. Though 
an experimental group reports a pygmy resonance at 3 
— 4 MeV 3l| , we see no indication of one in our calcula- 
tion. It has been suggested that in spherical Sn isotopes 
such resonances involve configurations beyond the natu- 
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FIG. 6: (Color online) Predicted IV E\ strength function in 
^^^Yb, with a folding width of 0.5 MeV (0.1 MeV for dis- 
crete states). 0~ and 1~ indicate K'^ components (the curve 
labeled K'^ = 1" includes the contribution of K'" = ~1~). 

ral ambit of the QRPA but the issue is unresolved, 

and in deformed nuclei has not been systematically in- 
vestigated. 

In summary, we have developed a box-based Skyrme- 
QRPA code for axially-symmetric even-even nuclei aimed 
at calculations throughout the isotopic chart. We rigor- 
ously tested the accuracy of the code in Mg isotopes, 



paying particular attention to spurious states. Though 
the code does not do quite so well with spurious states 
in heavier nuclei under the restrictions on space size that 
we currently impose, the agreement of energy weighted 
sums with sum rules indicates that a) in channels with 
no spurious modes (e.g. K'" = 2+) our code is quite ac- 
curate, b) in the K'^ = 0^ channel, spurious admixtures 
are negligible, and c) in the IS A'^^^O", and 1~ chan- 
nels (and the 1+ channel in the future), admixtures are 
larger but can be effectively removed. In addition, calcu- 
lations with larger spaces should be possible. Our imme- 
diate plans are to systematically investigate the ability 
of modern density functionals, together with the QRPA, 
to describe /3 and 7 vibrations in the rare-earth nuclei. 
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